Long Live the Disk: Lifetimes of Protoplanetary Disks in Hierarchical Triple-star Systems and a Possible Explanation for HD 98800 B

Our subject of study is a disk in a hierarchical triple-star system, that is, a circumbinary disk affected by an external stellar companion. Figure 1 shows a schematic view of this configuration that, due to the axisymmetrical nature of our model, will always consider circular and coplanar orbits for the inner binary, the circumbinary disk, and the outer stellar companion. Moreover, despite the disk being subject to a time-dependent gravitational potential due to the stellar orbital motions, we consider the disk to be rotating in a gravitational potential given by the total mass of the inner binary system and neglecting the gravitational perturbations of the external companion. We will discuss the limitations of these considerations in Section 6.

Zoom In Zoom Out Reset image size

To compute the vertical structure of a disk in this environment, we assume an axisymmetric, irradiated disk in hydrostatic equilibrium. We follow Guilera et al. (2017, 2019, 2020), who used the same methodology described in Alibert et al. (2005) and Migaszewski (2015) for a circumstellar case, but considering extra tidal heating terms in the energy equation due to the effects of the inner binary and the external companion (Lodato et al. 2009; Fontecilla et al. 2019). The complete set of equations is given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial z}=-\rho {{\rm{\Omega }}}^{2}z,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial F}{\partial z}=\displaystyle \frac{9}{4}\rho \nu {{\rm{\Omega }}}^{2}+{D}_{{{\rm{\Lambda }}}_{{\rm{I}}}}+{D}_{{{\rm{\Lambda }}}_{\mathrm{II}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial z}={\rm{\nabla }}\displaystyle \frac{T}{P}\displaystyle \frac{\partial P}{\partial z},\end{eqnarray} \tag{ 3 }$$

where P, ρ, F, T, and z represent the pressure, density, radiative heat flux, temperature and vertical coordinate of the disk, respectively. Here $\nu =\alpha {c}_{s}^{2}/{\rm{\Omega }}$ is the viscosity (Shakura & Sunyaev 1973), where α is a dimensionless parameter, ${c}_{s}^{2}=P/\rho$ is the square of the locally isothermal sound speed, ${\rm{\Omega }}=\sqrt{{{GM}}_{{\rm{I}}}/{R}^{3}}$ is the Keplerian frequency at a given radial distance R from the central binary, and M_I = M₁ + M₂ is its mass. We consider that the heat in the circumbinary disk is vertically transported by radiation or convection following the standard Schwarzschild criterion, as described in Guilera et al. (2019).

The terms ${D}_{{{\rm{\Lambda }}}_{{\rm{I}}}}=({{\rm{\Omega }}}_{{\rm{I}}}-{\rm{\Omega }}){{\rm{\Lambda }}}_{{\rm{I}}}\rho$ and ${D}_{{{\rm{\Lambda }}}_{\mathrm{II}}}=({{\rm{\Omega }}}_{\mathrm{II}}-{\rm{\Omega }}){{\rm{\Lambda }}}_{\mathrm{II}}\rho$ in Equation (2) represent the tidal heating dissipated in the form of radiation due to the inner binary and the outer stellar companion, respectively. Here Ω_I is the Keplerian frequency at the inner binary separation, a_I, given by ${{\rm{\Omega }}}_{{\rm{I}}}=\sqrt{G({M}_{{\rm{I}}})/{a}_{{\rm{I}}}^{3}}$, and Ω_II is the Keplerian frequency at the separation between the inner binary system and the external companion star, a_II, computed as ${{\rm{\Omega }}}_{\mathrm{II}}=\sqrt{G({M}_{\mathrm{II}})/{a}_{\mathrm{II}}^{3}}$, where M_II = (M₁ + M₂) + M₃, the total mass of the stellar system.

The quantities Λ_I and Λ_II represent the torques due to the inner binary and outer stellar companion, respectively. The classical expression for Λ proposed by Armitage & Natarajan (2002) ⁹ for q ≪ 1 could be overestimating the tidal torque contribution on the disks outside the Lindblad resonances, which may modify the surface density profiles (Tazzari & Lodato 2015). Thus, following Tazzari & Lodato (2015) and Fontecilla et al. (2019), we consider an exponential cutoff for the tidal torques outside this region, computed as

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{{\rm{I}}}(R)=\displaystyle \frac{f}{2}{q}_{{\rm{I}}}^{2}{{\rm{\Omega }}}^{2}{R}^{2}{\left(\displaystyle \frac{{a}_{{\rm{I}}}}{{{\rm{\Delta }}}_{{\rm{I}}}}\right)}^{4}{e}^{-{\left(\displaystyle \frac{R-{R}_{\mathrm{OM}}}{{W}_{\mathrm{OM}}}\right)}^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\mathrm{II}}(R)=-\displaystyle \frac{f}{2}{q}_{\mathrm{II}}^{2}{{\rm{\Omega }}}^{2}{R}^{2}{\left(\displaystyle \frac{R}{{{\rm{\Delta }}}_{\mathrm{II}}}\right)}^{4}{e}^{-{\left(\displaystyle \frac{R-{R}_{\mathrm{IM}}}{{W}_{\mathrm{IM}}}\right)}^{2}},\end{eqnarray} \tag{ 5 }$$

where f is a dimensionless normalization parameter that can adopt different values between 0.001 and 1 (Armitage & Natarajan 2002; Alexander 2012; Vartanyan et al. 2016; Shadmehri et al. 2018). The variable q_I represents the mass ratio between M₂ and M₁, and q_II is the mass ratio between M₃ and (M₁ + M₂). Finally, ${{\rm{\Delta }}}_{{\rm{I}}}=\max ({R}_{\mathrm{Hill}}^{{\rm{I}}},{{\rm{H}}}_{{\rm{g}}},| R-{a}_{{\rm{I}}}| )$ and ${{\rm{\Delta }}}_{\mathrm{II}}=\max ({R}_{\mathrm{Hill}}^{\mathrm{II}},{{\rm{H}}}_{{\rm{g}}},| R-{a}_{\mathrm{II}}| )$, where ${R}_{\mathrm{Hill}}^{{\rm{I}}}={a}_{{\rm{I}}}{({q}_{{\rm{I}}}/3)}^{1/3}$ is the Hill radius of the secondary star in the binary system, ${R}_{\mathrm{Hill}}^{\mathrm{II}}={a}_{\mathrm{II}}{({q}_{\mathrm{II}}/3)}^{1/3}$ is the Hill radius of the third star, and H_g is the height scale of the disk. As in Fontecilla et al. (2019), we consider R_OM = 1.59a_I and ${R}_{\mathrm{IM}}=0.63{a}_{\mathrm{II}}$ to be the radii of the outermost and innermost Lindblad resonances and W_OM =75H_g and ${W}_{\mathrm{IM}}=370{{\rm{H}}}_{{\rm{g}}}$ to be the widths of the Gaussian smoothing.

Note that, if the case under study is that of a circumstellar disk, Λ_I = Λ_II = 0, while if it is a circumbinary disk, only Λ_II = 0, and if it is a circumstellar disk affected by an external star companion, only Λ_I = 0.

The boundary conditions at the surface of the disk H, P_s = P(z = H), F_s = F(z = H), and T_s = T(z = H), considering the effects of the inner binary and outer star, are given by

$$\begin{eqnarray}&&{P}_{s}=\displaystyle \frac{{{\rm{\Omega }}}^{2}H{\tau }_{\mathrm{ab}}}{{\kappa }_{s}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{l}{F}_{s}=\displaystyle \frac{3{\dot{M}}_{\mathrm{st}}{{\rm{\Omega }}}^{2}}{8\pi }\\ \qquad +\displaystyle \frac{(({{\rm{\Omega }}}_{{\rm{I}}}-{\rm{\Omega }}){{\rm{\Lambda }}}_{{\rm{I}}}+({{\rm{\Omega }}}_{\mathrm{II}}-{\rm{\Omega }}){{\rm{\Lambda }}}_{\mathrm{II}}){\dot{M}}_{\mathrm{st}}{\rm{\Omega }}\mu {m}_{{\rm{H}}}}{6\pi \alpha {k}_{{\rm{B}}}{\left({T}_{\mathrm{irr}}^{4}+{T}_{s}^{4}\right)}^{1/4}},\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}\begin{array}{rcl}0 & = & 2\sigma \left({T}_{s}^{4}+{T}_{\mathrm{irr}}^{4}-{T}_{{\rm{b}}}^{4}\right)-\displaystyle \frac{9\alpha {\rm{\Omega }}k{({T}_{s}^{4}+{T}_{\mathrm{irr}}^{4})}^{1/4}}{8{\kappa }_{s}\mu {m}_{{\rm{H}}}}\\ & & -{F}_{s}-({{\rm{\Omega }}}_{{\rm{I}}}-{\rm{\Omega }})\displaystyle \frac{{{\rm{\Lambda }}}_{{\rm{I}}}}{2{\kappa }_{s}}-({{\rm{\Omega }}}_{\mathrm{II}}-{\rm{\Omega }})\displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{II}}}{2{\kappa }_{s}},\end{array}\end{eqnarray} \tag{ 8 }$$

where τ_ab = 0.01 is the optical depth, T_b = 10 K is the background temperature, and ${\dot{M}}_{\mathrm{st}}$ is the equilibrium accretion rate. In addition to the surface boundary conditions, due to symmetry, we consider that the heat flux in the midplane of the disk should be F(z = 0) = F₀ = 0.

Since we now have an inner binary system, the temperature associated with the stellar irradiation is considered as

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}^{4}={T}_{{\mathrm{irr}}_{1}}^{4}+{T}_{{\mathrm{irr}}_{2}}^{4};\end{eqnarray} \tag{ 9 }$$

i.e., we are adding up the fluxes due to each star, making the approximation that both are located at their center of mass, where

$$\begin{eqnarray}\begin{array}{rcl}{T}_{{\mathrm{irr}}_{i}} & = & {T}_{i}\left[\displaystyle \frac{2}{3\pi }{\left(\displaystyle \frac{{R}_{{\star }_{i}}}{R}\right)}^{3}+\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{R}_{{\star }_{i}}}{R}\right)}^{2}\right.\\ & & {\left.\times \left(\displaystyle \frac{H}{R}\right)\left(\displaystyle \frac{d\mathrm{log}H}{d\mathrm{log}R}-1\right)\right]}^{0.5},{\rm{i}}=1,2,\end{array}\end{eqnarray} \tag{ 10 }$$

and ${R}_{{\star }_{i}}$ and T_i are the radius and effective temperature of the primary (i = 1) and secondary (i = 2) stars of the inner binary system, R is the radial coordinate, and $d\mathrm{log}H/d\mathrm{log}R=9/7$. The individual stellar masses, with their corresponding stellar radii and temperatures, are obtained from Baraffe et al. (2015). Note here that our model does not consider the possible effect of self-shadowing if the disk becomes nonflared or the possible irradiation from the third star.

Equations (1)–(3) of the vertical structure are computed using a shooting method with a Runge–Kutta–Fehlberg integrator, while Equations (6)–(8) are solved by a multidimensional Newton–Raphson algorithm. Solving these equations allows us to obtain the mean viscosity $\overline{\nu }({{\rm{\Sigma }}}_{{\rm{g}}},R)$, where Σ_g is the gas surface density, required to solve the radial evolution of the disk.

It is worth noting that the detailed computation of the vertical structure and the thermodynamical quantities of circumbinary disks in hierarchical triple-star systems can be particularly useful for future transfer radiative calculations.

In order to compute the time evolution of the gas surface density, we solve the classical 1D diffusion equation (Pringle 1981), including the terms corresponding to the torques injected by the binaries, already described in the previous subsection,

$$\begin{eqnarray}\begin{array}{rcl}\frac{{\rm{\partial }}{{\rm{\Sigma }}}_{{\rm{g}}}}{{\rm{\partial }}t} & = & \frac{3}{R}\frac{{\rm{\partial }}}{{\rm{\partial }}R}\left[{R}^{1/2}\frac{{\rm{\partial }}}{{\rm{\partial }}R}(\bar{\nu }{{\rm{\Sigma }}}_{{\rm{g}}}{R}^{1/2})-\frac{2{{\rm{\Sigma }}}_{{\rm{g}}}({{\rm{\Lambda }}}_{{\rm{I}}}+{{\rm{\Lambda }}}_{\mathrm{II}})}{3{\rm{\Omega }}}\right]\\ & & +{\dot{{\rm{\Sigma }}}}_{{\rm{w}}}(R),\end{array}\end{eqnarray} \tag{ 11 }$$

where t is the temporal coordinate, $\overline{\nu }$ is the mean viscosity at the midplane of the disk, Σ_g is the gas surface density, and ${\dot{{\rm{\Sigma }}}}_{{\rm{w}}}(R)$ is the sink term due to the X-ray photoevaporation. The latter is computed following the analytical prescriptions derived by Owen et al. (2012) from radiation-hydrodynamic models (Owen et al. 2011, 2012). The total mass-loss rate depends on the X-ray luminosity of the central star, which in our case is actually a binary. Therefore, we use

$$\begin{eqnarray}&&{L}_{{\rm{X}}}={L}_{{\rm{X}}}^{1}+{L}_{{\rm{X}}}^{2},\end{eqnarray} \tag{ 12 }$$

where ${L}_{{\rm{X}}}^{1}$ and ${L}_{{\rm{X}}}^{2}$ are the X-ray luminosities of the stars of the inner binary system, which can be computed following Preibisch et al. (2005),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({L}_{{\rm{X}}}^{i}[\mathrm{erg}\,{{\rm{s}}}^{-1}]) & = & 30.37\\ +1.44\mathrm{log}\left(\displaystyle \frac{{M}_{i}}{\,{M}_{\odot }}\right),{\rm{i}} & = & 1,2,\end{array}\end{eqnarray} \tag{ 13 }$$

who derived this correlation between X-ray luminosity and stellar mass for young ($5\lt \mathrm{log}\tau [\mathrm{yr}]\lt 9.5$) low-mass stars (M_⋆ ≤ 2 M_⊙). Then, as described in Appendix B of Owen et al. (2012), for primordial disks (or disks without holes), the total mass-loss rate is given by

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{w}}} & = & 6.25\times {10}^{-9}{\left(\displaystyle \frac{{M}_{{\rm{I}}}}{1\,{M}_{\odot }}\right)}^{-0.068}\\ & & \times {\left(\displaystyle \frac{{L}_{{\rm{X}}}}{{10}^{30}\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{1.14}\,\,{M}_{\odot }{\mathrm{yr}}^{-1},\end{array}\end{eqnarray} \tag{ 14 }$$

while for disks with holes or transition disks, it is given by

$$\begin{eqnarray}\begin{array}{rcl}{\dot{M}}_{{\rm{w}}} & = & 4.8\times {10}^{-9}{\left(\displaystyle \frac{{M}_{{\rm{I}}}}{1\,{M}_{\odot }}\right)}^{-0.148}\\ & & \times {\left(\displaystyle \frac{{L}_{{\rm{X}}}}{{10}^{30}\mathrm{erg}\,{{\rm{s}}}^{-1}}\right)}^{1.14}\,\,{M}_{\odot }{\mathrm{yr}}^{-1}.\end{array}\end{eqnarray} \tag{ 15 }$$

These approximations were derived for circumstellar disks that can present both kinds of morphologies at different epochs. However, circumbinary disks always present inner cavities; thus, we only compute the X-ray mass-loss rate following Equation (15).

From the normalized radial mass-loss profiles derived by Owen et al. (2012) and presented in their Appendix B, and by computing ${\dot{M}}_{{\rm{w}}}$ (from Equation (15) in our case), we can compute ${\dot{{\rm{\Sigma }}}}_{{\rm{w}}}(R)$ that verifies

$$\begin{eqnarray}&&{\dot{M}}_{{\rm{w}}}=\int 2\pi r{\dot{{\rm{\Sigma }}}}_{{\rm{w}}}(R){dR}.\end{eqnarray} \tag{ 16 }$$

As in Rosotti & Clarke (2018), we are not considering the possible contribution from the external companion to the photoevaporation rate of the circumbinary disk. The choice to only consider X-ray photoevaporation in our model is due to the fact that this type of irradiation has yielded much higher photoevaporation rates than the extreme-UV (EUV) rates (Ercolano & Owen 2010; Owen et al. 2011; Picogna et al. 2019). Moreover, Kunitomo et al. (2021) recently showed that this is the predominant disk dispersal mechanism for low-mass stars (<2.5 M_⊙), compared to EUV and far-UV photoevaporation (see their Figure 11).

An important point to mention here is that, in the 1D approach, the form of Equations (4) and (5) completely prevents viscous gas accretion onto the stars. However, as shown by Artymowicz & Lubow (1996) and many others (e.g., Cuadra et al. 2009; Dunhill et al. 2015; Tang et al. 2017), accretion onto the stars still happens via tidal streams, a mechanism that can only be captured in 2D and 3D simulations. This mass-loss rate is not constant but modulated by the binary orbit. While MacFadyen & Milosavljević (2008) estimated that this accretion could be an average of ∼10% of the one expected in the absence of the binary, D'Orazio et al. (2013) and Farris et al. (2014) suggested that it could reach up to 30%–60%. The results of Ragusa et al. (2016) are in agreement with these previous works when the disk is thick. However, these authors found that it is possible to suppress all mass accretion if the disk becomes sufficiently thin. Due to the uncertainties on the amount of this mass loss, our model allows that a certain fraction of the expected accretion rate is accreted by the inner binary. In a similar way as Alexander (2012), we compute the accretion rate ${\dot{M}}_{{\rm{B}}}$ as ${\dot{M}}_{{\rm{B}}}=3\pi \nu {{\rm{\Sigma }}}_{{\rm{g}}}$ at every radial bin between R_c and 2R_c (where R_c is the radius of the inner cavity) and adopt the maximum value of ${\dot{M}}_{{\rm{B}}}$ in that region, considering it as an upper limit.

The disk radial structure is solved using an implicit Crank–Nicholson method from a_I to a_II using 2000 radial bins logarithmically equally spaced. We set boundary conditions of zero torque, which is equivalent to having zero density at the boundaries. We note here that the boundary conditions are applied to the radial bins corresponding to a_I and a_II. However, as we show in the next sections, the tidal torques included in the model naturally truncate the disk at different locations. A validation of our 1D+1D model for the case of circumbinary disks is presented in Appendix A, comparing our results to those of Vartanyan et al. (2016).

It has been previously shown that circumbinary disks dissipate at a lower rate than circumstellar disks (Alexander 2012). However, what happens to disks in triple-star systems? To address this question, we performed two simulations: our fiducial case T1 (see Table 1), which represents a circumbinary disk in a triple hierarchical system, and B1, which is exactly the same but ignoring the external star, thus representing an isolated circumbinary disk. We recall that, in terms of our model, the latter is achieved by setting Λ_II = 0.

4.1.1. Without Accretion by Tidal Streams

In order to focus on how the disk density profile changes with time, we first do not allow any gas accretion via streams to occur in these simulations. Recall from the end of Section 2.2 that, in 1D models, the form of the torque given by Equation (4) prevents any mass loss by viscous accretion unless accretion by tidal streams is allowed. Thus, the disk only loses mass by photoevaporation. Since the mass-loss rate due to the X-ray photoevaporative winds only depends on the X-ray luminosity of the central object (see Equation (15)), and since T1 and B1 systems have the same disk mass and inner stellar pair, the dissipation timescale for both disks is exactly the same. However, the shapes of the density profiles are quite different. In terms of our model, a gas disk has completely dissipated when its mass in our simulations is lower than 1 × 10⁻⁶ M_⊙. If, as in this case, the disk mass only evolves due to X-ray photoevaporation, the dissipation timescale can be directly estimated without the need of performing simulations by computing L_X in Equation (12) and then ${\dot{M}}_{{\rm{w}}}$ in Equation (15). For both B1 and T1, q_I = 1 and M₁ = M₂ = 0.5 M_⊙; thus, we obtain ${\dot{M}}_{{\rm{w}}}\approx 8.9\times {10}^{-9}\,\,{M}_{\odot }$ yr⁻¹, and both disks dissipate in τ ≈ 5.6 Myr.

Figure 2 shows the time evolution of the gas surface density (top), temperature (middle), and aspect ratio (bottom) profiles for the disk of our fiducial case (T1) in a hierarchical triple-star system (left column) compared to the same disk in a binary system (B1; right column) and, for both cases, without any mass loss by viscous accretion.

Zoom In Zoom Out Reset image size

The differences between the disk profiles in the triple and binary systems are remarkable. At early times (0.5 Myr of evolution; red lines) and due to the injection of angular momentum by the inner binary, the inner radius of the gas disk expands, creating an inner cavity of up to ∼1.65 au (∼3.30a_I, as is expected from previous studies; e.g., Artymowicz & Lubow 1994) for both cases, while the outer radius presents different values. In T1, the effect of the outer stellar companion located at 100 au prevents disk expansion beyond ∼35 au, which represents ∼0.35a_II and is in agreement with previous studies (Papaloizou & Pringle 1977; Pichardo et al. 2005). This situation is different for the circumbinary case, B1, where the disk can freely expand beyond ∼100 au. As a natural consequence of this expansion, the maximum value reached by the gas surface density profile is lower than that for the one in the triple-star system by about a factor of ∼0.5, and this difference increases with time, although the mass for both cases evolves at the same rate (i.e., each comparable, same-color profile between T1 and B1 presents the same amount of mass).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the gas density profiles in T1 evolve slowly in a confined region without significantly changing their shape, in B1, the gas disk is pushed further out, generating a growing cavity as a simultaneous consequence of photoevaporation and disk expansion due to angular momentum conservation (see Figure 3 for a comparison between the evolution of the inner cavity sizes). In addition, the gas surface density profiles in B1 decay faster than in T1, thus decreasing the local available mass of gas inside ∼35 au. As an example, at 4 Myr, the gas surface density at 10 au is ∼10 g cm⁻² for B1, while it is ∼200 g cm⁻² for T1. This situation could favor gas-giant planet formation in triple-star systems, in which a protoplanet can spend more time embedded in the gas disk compared to a circumbinary disk. However, it is important to note that higher gas surface densities can also enhance inward migration of low- and intermediate-mass planets.

Zoom In Zoom Out Reset image size

An interesting point here is that, as a consequence of the higher gas surface density profiles for T1, the midplane temperature and aspect ratio profiles of the same disk are higher than those of B1. These quantities could play an important role in the solid-component evolution of these protoplanetary disks and thus in the process of planet formation, since higher gas densities also imply higher dust densities. As an example, the location of the water-ice line (∼170 K), beyond which an increment in the solid surface density is expected, is always further out in T1 compared to B1, and it also lasts longer, as can be seen in Figure 3. This situation may lead to the formation of ice cores in T1 for longer timescales and at larger distances than in B1. In addition, higher aspect ratios reduce the migration rate of low- and intermediate-mass planets (e.g., Paardekooper et al. 2011; Jiménez & Masset 2017) and increase the pebble isolation mass (e.g., Lambrechts et al. 2014). These two effects could help planet formation in hierarchical triple-star systems more than in binaries.

On the other hand, we note that while the temperature radial profiles are similar in both cases, i.e., they decrease as we move away from the central binary, the aspect ratio profiles present more differences. This is particularly important for self-shadowing, although we reiterate here that this effect is not included in the disk thermodynamics of our model. In any case, self-shadowing would be stronger in the case of the triple than the circumbinary-star system. In the latter case, the decrease in the aspect ratio could generate an effect of self-shadowing until ∼10 au, leading to a cooler disk between ∼2 and ∼10 au than previously computed. This effect would be negligible after 2.5 Myr. For the disk in the triple-star system, self-shadowing would affect it until ∼20 au at early times and persist until ∼4.5 Myr of the disk evolution.

4.1.2. With Accretion by Tidal Streams

We now describe the global results of the simulations of Table 1 for disks that evolve due to photoevaporation and viscous accretion by tidal streams with an efficiency of 30%.

4.2.1. Dependence on the Stellar Mass Ratios

4.2.2. Dependence on the Semimajor Axes

As mentioned above, HD 98800 is a quadruple-star system. However, in order to model the circumbinary disk around Ba–Bb using PLANETALP-B, we consider the outer binary Ab–Aa as a single body. This is a good approximation of its gravitational effect given the large separation between the binaries, and also the orbit of this component has not been resolved yet. Moreover, due to the axial symmetry of our model, we need to adopt an important simplification: circular and coplanar orbits for both star systems and the disk. We discuss these simplifications in Section 6.

We model the time evolution of the disk around HD 98800 B by considering stars of masses 0.7 and 0.6 M_⊙ for the system Ba–Bb, so their mass ratio is q_Ba–Bb = 0.857, with a total mass of 1.3 M_⊙ for the system Aa–Ab. We consider two different setups for the orbital parameters. Setup 1 adopts the observed semimajor axes of the system, namely, a_Ba–Bb = 0.5 and a_A−B = 45 au, but in a circular and coplanar configuration due to the geometrical limitations of our model. For setup 2, we adopt the extreme values a_Ba–Bb = 0.9 and a_A−B = 27 au, which correspond to the maximum and minimum expansion of the Ba–Bb and B–A orbits (apocenter and pericenter), considering their observed eccentricities of 0.78 and 0.4, respectively, but again, in a circular and coplanar configuration. This range of parameters should capture some of the effects we miss by not having eccentric orbits in the models, allowing us to get realistic estimates of the minimum dissipation timescale of the system.

Unlike in Section 4, we now fix the stellar and orbital parameters of the system and consider different disk parameters, such as different values for M_d and α and different accretion efficiency fractions , shown in Table 2. As in our previous simulations of Section 4, we keep a fixed exponent for the initial density profile of γ = 1 and a fixed characteristic radius of R₀ = 39 au.

Kastner et al. (2004) were able to measure the X-ray luminosity of HD 98800 B, which is ${L}_{{\rm{X}}}^{\mathrm{Ba}\mbox{--}\mathrm{Bb}}\approx 1.4\times {10}^{29}$ erg s⁻¹. This value is an order of magnitude lower than the one we would get with our approach, ${L}_{{\rm{X}}}^{\mathrm{Ba}\mbox{--}\mathrm{Bb}}={L}_{{\rm{X}}}^{\mathrm{Ba}}+{L}_{{\rm{X}}}^{\mathrm{Bb}}\approx 2.53\,\times {10}^{30}$ erg s⁻¹ following Equations (12) and (13), which has important consequences for the disk evolution. As explained in Section 4.1, if only X-ray photoevaporation is considered, and no mass loss due to accretion by tidal streams is allowed, the dissipation timescale can be directly estimated without the need of performing simulations. Using ${L}_{{\rm{X}}}^{\mathrm{Ba}\mbox{--}\mathrm{Bb}}$ as measured by Kastner et al. (2004), we obtain disk lifetimes much greater than 10 Myr for disk masses of M_d = 0.1, 0.05, and 0.01 M_⊙, while, using L_X as computed with Equations (12) and (13), we obtain disk lifetimes of 7.5, 3.7, and 0.7 Myr for the same disk masses. The latter values can be considered as lower limits for the HD 98800 B disk lifetime.

Taking advantage of the fact that this system has a measured X-ray luminosity, we then run the simulations only considering ${L}_{{\rm{X}}}^{\mathrm{Ba}\mbox{--}\mathrm{Bb}}$ as estimated by Kastner et al. (2004). Due to the uncertainties in the mass-loss rates due to stream accretion, we also consider different values for the accretion efficiency , as shown in Table 2. In Appendix B, we analyze the case in which the accretion efficiency fraction is not set as a constant value but rather as a function of the aspect ratio of the disk, as proposed by Ragusa et al. (2016).

Figure 6 shows the time evolution of the gas surface density (left column), temperature (middle column), and aspect ratio (right column) profiles every 0.1 Myr, represented by the color scale, for sims 4 (top), 5 (middle), and 6 (bottom) of setup 1 (see Table 2) with an accretion efficiency of 30% ( = 0.3). These are the results for a disk of M_d = 0.05 M_⊙ with different α values in the rows. As mentioned before, our simulations end when the disk has completely dissipated, defined as a gas disk mass lower than 10⁻⁶ M_⊙, or when it reaches 10 Myr.

Zoom In Zoom Out Reset image size

The viscosity parameter α, which is the only difference between the simulations of Figure 6, plays an important role in two different aspects. First, it is responsible for shaping the gas surface density, temperature, and disk aspect ratio profiles and, consequently, determining the size and width of the disk. Second, it directly affects the dissipation timescales when some accretion is allowed via streams.

Regarding sizes, the higher the value of α, the wider the disk, which expands in both directions, toward the central binary and outer star. The boundaries for each case, R_c and R_t, can be appreciated in Table 2. These values correspond to the gas disk profile at ∼0.05 Myr of evolution and measured at Σ_g = 10⁻³ g cm⁻² and were computed for the case in which no viscous accretion via streams is allowed. However, they do not change significantly for the cases in which it is allowed; thus, for simplicity, we only show in Table 2 the sizes of the disks for the = 0 efficiency cases for setups 1 and 2.

The time evolution of the gas surface density profiles for disks with M_d = 0.01 and 0.1 M_⊙ behave in a similar way to the ones described before for M_d = 0.05 M_⊙ but with different dissipation timescales.

Table 2 also shows the dissipation timescales for each of the performed simulations for both setup 1 and setup 2. Those simulations in which the disk dissipates only after the minimum estimated age of HD 98800, 7 Myr, are highlighted in bold for setups 1 and 2. However, only those additionally marked with an asterisk also fit the estimated mass of ∼5 M_J of the disk around HD 98800 B (Ribas et al. 2018). The latter can also be seen in Figure 7, which shows the time evolution of the mass of the disk for all setup 1 simulations. The bottom panel shows the results for sims 1, 2, and 3 with a disk of M_d = 0.01 M_⊙; the middle panel is for sims 4, 5, and 6 with a disk of M_d = 0.05 M_⊙; and the top panel is for sims 7, 8, and 9 with a disk of M_d = 0.1 M_⊙. Green curves represent simulations with α = 0.0001, purple with α = 0.001, and orange with α = 0.01. Solid lines represent simulations with accretion efficiencies of 0%, short-dashed lines are 10%, medium-dashed lines are 30%, and long-dashed lines are 50%. The horizontal gray dashed line marks the current estimated mass of HD 98800 B (Ribas et al. 2018), and the shaded area is its probable age. The curves marked with black arrows are the ones highlighted with an asterisk in Table 2.

Zoom In Zoom Out Reset image size

From Table 2 and Figure 7, it can be deduced that, despite the complex configuration of HD 98800, and assuming our simplified model for HD 98800, we found several cases that would explain the existence of the massive disk after ∼7 Myr, even with high mass accretion efficiencies via streams and relatively low initial masses. However, moderate-to-low viscosities are needed for the low initial mass cases. It is interesting to remark, though, that the low-viscosity values are the ones that better reproduce some of the ring structures in the DSHARP disks (Dullemond et al. 2018).

As shown by the results in Table 2, for both setup 1 and setup 2, highly viscous disks that are allowed to lose mass via tidal streams with ≥10% efficiency dissipate in less than 3.81 Myr, even for the most conservative simulations (the most massive disk in setup 1).

Given the similar total stellar masses for systems A and B, the lack of a disk in system A is intriguing. We do know, however, that the X-ray luminosity of HD 98800 A is about four times higher than that of HD 98800 B (Kastner et al. 2004), which implies a shorter lifetime for a disk around it. Here we will test this idea quantitatively with our model.

We will assume for simplicity that the binary and initial disk in system A have the same parameters as those in system B. It is important to remark that, because system Aa–Ab was classified as an SB1 instead of SB2 (Torres et al. 1995), its actual mass ratio should be lower than that of B. As shown in Section 4.2.1, a smaller value of q_I leads to a faster disk dispersal. Therefore, our choice of q_A = q_B, even if not very realistic, is conservative, as it makes the disk in the A system last for longer.

We compute the time evolution of disks around binary A for orbital setup 1, with ${L}_{{\rm{X}}}^{{\rm{A}}}=4{L}_{{\rm{X}}}^{{\rm{B}}}$, and using the parameters of Table 2 that resulted in dissipation timescales longer than 7 Myr for system B (highlighted in bold). Table 3 presents the dissipation timescales for these new simulations. Those in bold represent disks that dissipated in less than 7 Myr, and those additionally marked with an asterisk are the ones that, together with the results of Table 2 (and assuming that both disks, A and B, initially have the same mass), better fit the current mass and age of HD 98800 B and the current nonexistence of a disk around system A. For the disk around system A, moderate- to low-mass disks are needed to allow its dissipation timescale in less than 7 Myr, the minimum age estimated for the quadruple system.

One of the most important hypotheses in our model, a consequence of its 1D nature, is the axisymmetry of the disk. This characteristic does not properly reproduce the details of the disk structure, which is clearly nonaxisymmetric due to the non-Keplerian potential of the central binary. Even if this approximation might not be the most appropriate to describe the shape of the cavity and the inner regions of circumbinary disks (Dunhill et al. 2015; Miranda et al. 2017; Thun et al. 2017; Thun & Kley 2018; Poblete et al. 2019) heavily affected by the inner binary, it works well in the outer regions, where the binary potential can be approximated by the potential caused by a point mass. In addition, in the particular case of HD 98800 B, we remark that Kennedy et al. (2019) showed that the disk is largely axisymmetrical.

The second important assumption in our model is the circular and coplanar nature of the orbits of the inner binary and the circumbinary disk and external companion with respect to the inner binary. On the one hand, binaries on circular orbits produce the most eccentric cavities (Thun et al. 2017), a feature that we cannot reproduce with our 1D+1D model. Although the development of an eccentric cavity affects gas accretion, which results to be not constant and modulated by the inner binary orbit, we still take this accretion, on average, into account. For our modeling of the disk around HD 98800 B, though, this is not an issue, as the observed binary is not actually circular.

On the other hand, although a full coplanar orbital configuration in a hierarchical triple-star system could seem unlikely at first sight, a very recent analysis of the system TWA3 that hosts a circumbinary disk shows that this could be the case (Czekala et al. 2021). Moreover, we do not expect significant changes in our results for the disk around HD 98800 B, which in fact presents a polar configuration with respect to the inner binary system Ba–Bb (Kennedy et al. 2019), since Franchini et al. (2019) showed that for polar configurations, the main effect on the disk is a reduced inner cavity due to weaker inner torques. It is possible, though, that the irradiation and photoevaporation mechanisms affecting the circumbinary disk are different due to its polar configuration, since the solid angle intercepted by the disk is different. However, models capable of predicting the effects of irradiation and photoevaporation on disks from binary-star systems have not been developed yet, even less if the binary-star systems have a polar orbit with respect to their disks. Thus, we could not predict if this configuration would contribute to a faster disk dispersal or not.

Photoevaporation from the central star is perhaps one of the most important mechanisms affecting mass removal in protoplanetary disks and, as shown previously, is clearly important for circumbinary disks. Disks in hierarchical triple-star systems could have, in addition to the irradiation from the inner binary system, an extra source of X-ray radiation: the external companion. However, as described by Rosotti & Clarke (2018), the irradiation from the external star could be somehow important only if their separation is lower than the radius at which X-ray heating is effective (which is about 100 au; Owen et al. 2012). For the particular case of HD 98800 B, then, the irradiation from the external binary system Aa–Ab could somehow affect its evolution. However, 3D radiation-hydrodynamics simulations that have not yet been developed would be necessary to evaluate this effect.

Another mechanism that could affect disks in general, contributing to a faster disk dispersal, is external photoevaporation from close stellar encounters (Winter et al. 2018) or due to OB massive stars nearby (Clarke 2007; Anderson et al. 2013; Winter et al. 2020). Shadmehri et al. (2018) considered this latter effect in the evolution of circumbinary disks and found that it can significantly affect the disk masses and lifetimes depending on the intensity of the external UV radiation flux. We plan to include this effect in future works. However, it is important to remark that, at least for the Hya association, external photoevaporation does not seem to be relevant, since even with an estimated age of ∼10 Myr, it still presents disks with detected gas, such as HD 98800 B (Ribas et al. 2018) and TWA 3, where ¹²CO and ¹³CO were detected for the first time (Czekala et al. 2021).

One of the most important results of this work is that disks in hierarchical triple-star systems can survive for several megayears, even when allowing gas accretion by tidal streams. This situation, together with the fact that surface density values in triples are, in general, higher than the ones of their circumbinary counterparts (see Figure 2), could favor planet formation even more than in circumbinary disks. In addition, the location of the ice line in triple-star systems is always beyond the corresponding one in the circumbinary disk case, as we show in Section 4.1.1. Thus, for equal initial mass disks and dust-to-gas ratios in circumbinary and triple disks, higher surface densities and larger initial locations of the ice lines could favor the formation of ice cores at larger distances. This comparison between disks in binaries and triples is analog to the one between circumstellar and circumbinary disks, where the latter have higher gas surface densities. The fact that planet formation in triples, as in binaries, could be favored at larger distances is in line with the idea that planets in circumbinary disks form far away and then migrate inward, parking near the stability limit (Pierens & Nelson 2007; Thun & Kley 2018), as may have happened for the discovered Kepler circumbinary planets.

A possible problem that circumbinary disks in triple-star systems could have in achieving planet formation is their truncation by the external companion. Even in the case of having a disk in a binary and triple with the same initial mass of gas and solids, the truncation of the disk in the triple avoids the continuous supply of dust and pebbles from the outer region of the disk that could be quickly lost due to radial drift, unlike what happens in circumbinary disks not affected by an external companion. This situation was recently studied by Zagaria et al. (2021), who found that the dust in circumstellar disks around one of the stars in a binary-star system drifts faster than in disks not affected by an outer stellar companion.

Nevertheless, an interesting situation could occur for disks with large a_I (see panels T6 and T7 of Figure 5) and low viscosities (see bottom row of Figure 6). In these cases, the gas surface density shows quasi-Gaussian radial profiles with a well-defined maximum. As was recently shown by Guilera et al. (2020), a maximum in the gas surface density is related to a pressure maximum, and it represents a preferential location for the dust/pebble trap, planetesimal formation by streaming instability, and planet formation by hybrid accretion of pebbles and planetesimals. Moreover, the well-defined pressure maximum can act as a planet migration trap, thus favoring planet formation.